#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <fstream>
#include <math.h>
#include <cstdlib>
#include <unistd.h>
#include <vector>
#include <algorithm>


using namespace std;

#include "point.cpp"
#include "spotobins.cpp"

#include <xpa.h>

#define NXPA 10


class vShape{

public:
  vector <point> p;
  point bp;
  double deltaMAX;
  double alpha;
  double wid;
  vector <spotoBin> bin;
  double diamExpWeigth=0;
  double expAbove=1;

  vShape(point _bp, double _deltaMAX, double _alpha, double _wid)
  {
    setParam(_bp, _deltaMAX, _alpha,_wid);
  }

  void setParam(point _bp, double _deltaMAX, double _alpha, double _wid)
  {
    bp=_bp;
    alpha=_alpha;
    wid=_wid;
    deltaMAX=_deltaMAX;
  }

  void setDiamExpWeigth(double _diamExpWeigth)
  {
    diamExpWeigth=_diamExpWeigth;
  }
  
  void setExpAbove(double _expAbove)
  {
    expAbove=_expAbove;
  }

  vector <point> getPointsAbove()
  {
    double d=0;
    long N=p.size();
    double y0, yu, yd;
    double numerator, denominator;

    vector <point> pp;
    
    for (long i=0; i<N; i++)
      {
	y0=alpha*(p[i].x-bp.x)+bp.y;
	yd=alpha*(p[i].x-bp.x+wid)+bp.y;
	yu=alpha*(p[i].x-bp.x-wid)+bp.y;
	if ((p[i].y>=y0) && (p[i].y<yu))
	  pp.push_back(p[i]);
      }

    return pp;
      
  }
  
  vector <point> getPointsBelow()
  {
    double d=0;
    long N=p.size();
    double y0, yu, yd;
    double numerator, denominator;

    vector <point> pp;
    
    for (long i=0; i<N; i++)
      {
	y0=alpha*(p[i].x-bp.x)+bp.y;
	yd=alpha*(p[i].x-bp.x+wid)+bp.y;
	yu=alpha*(p[i].x-bp.x-wid)+bp.y;
	if ((p[i].y<y0) && (p[i].y>=yd))
	  pp.push_back(p[i]);
      }

    return pp;
      
  }

  double getValueLeft()
  {
    double d=0;
    long N=p.size();
    double y0, yu, yd;
    double above=0, below=0;

    
    for (long i=0; i<N; i++)
      {
	y0=-alpha*(p[i].x-bp.x)+bp.y;
	yd=-alpha*(p[i].x-bp.x+wid)+bp.y;
	yu=-alpha*(p[i].x-bp.x-wid)+bp.y;
	if ((p[i].y>=y0) && (p[i].y<=yu)) // above
	  //	  above+=1;
	  above+=pow(p[i].y,diamExpWeigth);
	if ((p[i].y<y0) && (p[i].y>=yd)) // below
	  //	  below+=1;
	  below+=pow(p[i].y,diamExpWeigth);
      }

    double pp;
    if (below==0)
      return 0;
    else 
      return pow(above,expAbove)/(below);
  }
  
  double getValueRight()
  {
    double d=0;
    long N=p.size();
    double y0, yu, yd;
    double above=0, below=0;

    
    for (long i=0; i<N; i++)
      {
	y0=alpha*(p[i].x-bp.x)+bp.y;
	yd=alpha*(p[i].x-bp.x-wid)+bp.y;
	yu=alpha*(p[i].x-bp.x+wid)+bp.y;
	if ((p[i].y>=y0) && (p[i].y<=yu)) // above
	  //	  above+=1;
	  above+=pow(p[i].y,diamExpWeigth);
	if ((p[i].y<y0) && (p[i].y>=yd)) // below
	  //	  below+=1;
	  below+=pow(p[i].y,diamExpWeigth);
      }

    double pp;
    if (below==0)
      return 0;
    else 
      return pow(above,expAbove)/(below);
  }

  double getValueBoth() // this is not working yet !
  {
    double d=0;
    long N=p.size();
    double y0, yu, yd;
    double above=0, below=0;

    
    for (long i=0; i<N; i++)
      {
	y0=alpha*(p[i].x-bp.x)+bp.y;
	yd=alpha*(p[i].x-bp.x-wid)+bp.y;
	yu=alpha*(p[i].x-bp.x+wid)+bp.y;
	if ((p[i].y>=y0) && (p[i].y<=yu)) // above
	  //	  above+=1;
	  above+=pow(p[i].y,diamExpWeigth);
	if ((p[i].y<y0) && (p[i].y>=yd)) // below
	  //	  below+=1;
	  below+=pow(p[i].y,diamExpWeigth);
      }

    double pp;
    if (below==0)
      return 0;
    else 
      return pow(above,expAbove)/(below);
  }
  
  /* this is bin2ascii output on april 2017
000001                Ceres  3.34 0.12  2.7671  0.1162  0.1676  973.890  13.310 0.087 0.003   8.95 C-----------G---C------------- 000000  0.00  0.00

   */
  int importPoint(string line)
  {
    double x, y;
    point pl;

    if (line.length()>0)
      {
	pl.x=atof(line.substr(39,8).c_str());
	pl.y=1/atof(line.substr(63,8).c_str());     
	p.push_back(pl);
	return 1;
      }
    else
      return 0;
  }

  vector <point> projPoints()
  {
    vector <point> pp;
    long N=p.size();
    
    for (long i=0; i<N; i++)
	pp.push_back(p[i].project(bp,alpha));

    return pp;
  }

  vector <point> projPointsAbove()
  {
    vector <point> pp=getPointsAbove();
    vector <point> ppp;
    
    long N=pp.size();
    
    for (long i=0; i<N; i++)
	ppp.push_back(pp[i].project(bp,alpha));

    return ppp;
  }

  point centerOfMass(vector <point> p)
  {
    long N=p.size();
    point c(0,0);
    
    for (long i=0; i<N; i++)
      {
	c.x+=p[i].x;
	c.y+=p[i].y;
      }
    
    c.x/=N;
    c.y/=N;
    
    return c;
  }

  double aveDista(vector <point> p)
  {
    long N=p.size();
    point c=centerOfMass(p);
    double d=0;

    for (long i=0; i<N; i++)
      d+=c.distance2(p[i]);
    
    d/=N;
    return sqrt(d);
  }
  /*
  void sortPoints()
  {
    sort(p.begin(), p.end(), [](const point* a, const point* b)
	 {
	   if(a < b) return -1;
	   if(a > b) return 1;
	   return 0;
	 });
  }
  */
  void sortPoints()
  {
    sort(p.begin(), p.end());
  }

  void makeBins(long N) // divide the range [bp.y, deltaMAX] in N bins
  {
    double deltaStep=(deltaMAX-bp.y)/N;
    spotoBin *pBin;

    for (long i=0; i<N; i++)
      {
	pBin = new spotoBin(bp.y+deltaStep*i, bp.y+deltaStep*(i+1), p);
	bin.push_back(*pBin);
	delete pBin;
      }
  }

  void sortBins()
  {
    sort(bin.begin(), bin.end());
  }
    
  void reBin(double std)
  {
    vector<spotoBin>::iterator b;
    vector <spotoBin> newBins;
    long NnewBins=0;
    
    do
      {
	for(vector<spotoBin>::iterator a = bin.begin(); a < bin.end()-1; a++)
	  {
	    b=a+1;
	    if ((b->getNpoints()-a->getNpoints())>std)
	      newBins.push_back(b->split());
	    if ((a->getNpoints()-b->getNpoints())>std)
	      newBins.push_back(a->split());	    
	  }
	bin.insert(bin.end(), newBins.begin(), newBins.end());
	NnewBins=newBins.size();
	newBins.clear();
	sortBins();
      } while (NnewBins>0);
  }

  void reBin()
  {
    double Nave=0;
    double Nstd=0;
    for(vector<spotoBin>::iterator a = bin.begin(); a < bin.end(); a++)
	Nave+=a->getNpoints();
    Nave/=bin.size();
    
    for(vector<spotoBin>::iterator a = bin.begin(); a < bin.end(); a++)
      Nstd=(Nave-a->getNpoints())*(Nave-a->getNpoints());
    Nstd/=bin.size();

    Nstd=sqrt(Nstd);

    reBin(Nstd);
  }
};


int main(int argc, char** argv)
{
  string line;
  //  point p(2.372,0);
  point p(2.372,0);
  vShape V(p,.15,-0.6,0.04);
  double wid=0.03;
  int right=0; // default is left side of the V
  int both=0; // count for a full V: i.e. left and right 
  
  opterr = 0;
  int c;
  int heatMap=1;
  int spotoBin=0;

  int Nbins=16;
  double a0=2.1;
  double a1=2.5;
  double K0=0.2;
  double K1=2.5;
  double aStep=0.01;
  double kStep=0.1;
  
  while ((c = getopt (argc, argv, "d:D:ba:A:k:K:e:w:p:rSMn:")) != -1)
    switch (c)
      {
      case 'K':
	K1=atof(optarg);
	break;
      case 'k':
	K0=atof(optarg);
	break;	
      case 'A':
	a1=atof(optarg);
	break;
      case 'a':
	a0=atof(optarg);
	break;
      case 'd':
	aStep=atof(optarg);
	break;
      case 'D':
	kStep=atof(optarg);
	break;
      case 'b':
	both=1;
        break;
      case 'r':
	right=1;
        break;
      case 'p':
	V.setDiamExpWeigth(atof(optarg));
        break;
      case 'e':
	V.setExpAbove(atof(optarg));
        break;
      case 'w':
	wid=atof(optarg);
        break;
      case 'S':
	spotoBin=1;
	heatMap=0;
        break;
      case 'M':
	spotoBin=0;
	heatMap=1;
        break;
      case 'n':
	Nbins=atoi(optarg);
	break;
      case '?':
        if (optopt == 'G')
          fprintf (stderr, "Option -%c requires an argument.\n", optopt);
        else if (optopt == 'p')
          fprintf (stderr, "Option -%c requires an argument.\n", optopt);
        else if (isprint (optopt))
          fprintf (stderr, "Unknown option `-%c'.\n", optopt);
        else
          fprintf (stderr,
                   "Unknown option character `\\x%x'.\n",
                   optopt);
        return 1;
      default:
        abort ();
      }

  
  while(cin)
    {
      getline(cin, line);
      if (V.importPoint(line)) {};
    }

  /*
  cout  << V.getValue() << endl;
  return 0;
  */

  /*  
  vector <point> pp=V.getPointsAbove();
  for (long i=0; i<pp.size(); i++)
    pp[i].print();

  cout << endl;
  cout << V.getValue() << endl;
  */
  
  //  vector <point> pp=V.projPointsAbove();
  // cout << V.aveDista(pp) << endl;

  if (spotoBin) {
    V.sortPoints();
    V.makeBins(Nbins);
    //  for (long i=0; i<V.bin.size(); i++) V.bin[i].print();
    
    //  V.reBin();
    //V.reBin();
    //V.reBin();
  
    //  for (long i=0; i<V.bin.size(); i++) V.bin[i].print();
    
    // return 0;
    for (long i=0; i<V.bin.size(); i++)
      {
	V.bin[i].getLeftMost().print();
      }
  
    return 0;
  }

  
  //THIS IS TO PLOT THE HEATMAP
  if (heatMap) {
    while ((a0+=aStep)<a1)
      {
	double alpha=K0;
	while((alpha+=kStep)<K1)
	  {
	    point p(a0,0);
	    V.setParam(p,0.5,alpha,wid/alpha);
	    if (both)
	      {
		cout << a0 << " " << alpha << " " << V.getValueBoth() << endl;
		continue;
	      }
	    if (right)	    
	      cout << a0 << " " << alpha << " " << V.getValueRight() << endl;
	    else
	      cout << a0 << " " << alpha << " " << V.getValueLeft() << endl;
	    //	  V.setParam(p,0.5,alpha,0.03/alpha);
	    //	  vector <point> pp=V.projPointsAbove();
	  }
	cout << endl;
      }
    return 1; // exit the heatMap mode
  }

  /*
  //THIS IS TO SEND THE HEATMAP TO DS9

  double minAlpha=0.1;
  double maxAlpha=2.5;
  double stepAlpha=0.02;
  long nAlpha=(maxAlpha-minAlpha)/stepAlpha;

  double minA0=2.1;
  double maxA0=2.5;
  double stepA0=0.002;
  long nA0=(maxA0-minA0)/stepA0;

  double* array = new double[nA0*nAlpha];
  
  double alpha=minAlpha;
  double a0=minA0;
  for (long ia0=0; ia0<nA0; ia0++)
    {
      a0=minA0+ia0*stepA0;
      point p(a0,0);
      for (long ialpha=0; ialpha<nAlpha; ialpha++)	  
	{
	  alpha=minAlpha+ialpha*stepAlpha;
	  if (right)
	    V.setParam(p,0.5,alpha,-wid/alpha);
	  else
	    V.setParam(p,0.5,-alpha,wid);	    
	  array[ia0+ialpha*nA0]=(double)V.getValue();
	  cout << a0 << " " << alpha << " " << array[ia0+ialpha*nA0] << " " << 2/alpha << endl;
	}
    }
  
  size_t  len=nA0*nAlpha*sizeof(array[0]);
  char *names[NXPA];
  char *messages[NXPA];
  int got;
  char param[256];
  //  sprintf(param,"array [xdim=%d,ydim=%d,bitpix=-16,endian=little]",nA0,nAlpha);
  sprintf(param,"array [xdim=%d,ydim=%d,bitpix=-64,endian=little]",nA0,nAlpha);
  std::cout << "PARAM TO XPA - ds9 " << param << std::endl;
  
  got = XPASet(NULL, "ds9", param, NULL, (char*)array, len, names, messages, NXPA);
  // error processing 
  for(int i=0; i<got; i++){
    if( messages[i] ){
      fprintf(stderr, "ERROR: %s (%s)\n", messages[i], names[i]);
    }
    if( names[i] )    free(names[i]);
    if( messages[i] ) free(messages[i]);
  }
  delete [] array;
  */
  
  return 0;     
}

/*
000084                 Klio    00084    9.32  0.15  2.3627  0.2362  0.1621  2.3622  0.1928  0.1637   79.000   4.867 0.053 0.017
000163              Erigone    00163    9.47  0.04  2.3673  0.1908  0.0839  2.3671  0.2096  0.0837   81.579   3.062 0.033 0.004
000249                 Ilse    00249   11.33  0.15  2.3781  0.2171  0.1671  2.3775  0.1789  0.1742   37.030   0.610 0.038 0.001
000282             Clorinde    00282   10.91  0.15  2.3392  0.0808  0.1570  2.3393  0.1026  0.1534   41.870   1.151 0.027 0.004
000284               Amalia    00284   10.05  0.11  2.3588  0.2215  0.1401  2.3584  0.1945  0.1546   57.840   0.720 0.051 0.001
000302             Clarissa    00302   10.89  0.15  2.4053  0.1118  0.0595  2.4057  0.1065  0.0584   34.440   0.480 0.068 0.002
000313             Chaldaea    00313    8.90  0.15  2.3755  0.1810  0.2020  2.3758  0.2325  0.2040   96.000   7.809 0.053 0.013
000336            Lacadiera    00336    9.76  0.15  2.2523  0.0955  0.0985  2.2518  0.0882  0.1071   69.000   3.364 0.046 0.005
000345            Tercidina    00345    8.71  0.10  2.3254  0.0612  0.1693  2.3253  0.0987  0.1799   99.000  11.469 0.059 0.012
*/
